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INTRODUCTION 


One  of  the  major  difficulties  in  the  formulation  of  effective  shell 
elements  has  been  identified  to  be  the  phenomenon  of  membrane  locking. 

Membrane  locking  occurs  in  curved  shell  elements  when  the  in-plane  displace¬ 
ment  approximation  is  not  of  higher  order  than  the  transverse  displacement 
approximation  and  inextensional  bending  of  the  shell  cannot  take  place. 
Inextentional  bending  is  an  important  mode  of  deformation,  and  when  an  element 
is  not  capable  of  representing  inextensional  bending,  parasitic  membrane 
energy  is  generated  in  many  modes  of  deformation.  In  the  same  manner  that 
parasitic  shear  causes  shear  locking,  this  spurious  membrane  energy  causes 
membrane  locking.  Membrane  locking  severely  reduces  the  rate  of  convergence 
of  shell  elements,  particularly  in  deep  shells  and  in  situations  where  the 
bending  of  the  shell  is  the  dominant  mode  of  deformation. 

In  this  report,  two  methods  for  eliminating  membrane  locking  in  curved 
shell  elements  are  presented.  The  first  method  is  a  strain  projection  method 
in  which  the  membrane  strains  are  corrected  so  that  inextensional  modes  of 
pure  bendiny  become  possible.  The  method  is  applied  both  to  a  curved  beam 
element  and  a  triangular  shell  element  in  which  the  flexural  behavior  is 
modeled  by  a  discrete  Kirchhoff  theory.  The  use  of  this  projection  method 
introduces  membrane-flexural  coupling  to  the  shell  element  and  modifies  the 
bendiny  stiffness  in  an  appropriate  fashion.  Results  have  been  obtained  with 
this  element  for  linear  analysis  of  static  response  of  deep  shell  structures 
and  for  nonlinear,  collapse  analysis  of  columns  and  cylindrical  panels.  The 
results  show  a  remarkably  rapid  rate  of  convergence. 


The  second  method  which  is  under  investigation  for  avoiding  membrane 
locking  is  the  use  of  uniform  reduced  quadrature  on  the  9-node  Lagrange 


element.  Results  obtained  in  the  previous  study  by  the  author  suggests  that 
if  2  x  2  quadrature  could  be  used  in  the  curved  shell  element,  membrane 
locking  as  well  as  shear  locking  could  be  avoided,  and  thus  very  rapid  rates 
of  convergence  would  be  achieved.  However,  the  use  of  uniform  reduced 
quadrature  is  accompanied  by  the  appearance  of  spurious  zero  energy  modes 
which  can  yield  meaningless  answers  to  solutions  for  certain  boundary 
conditions . 

In  this  investigation,  a  spurious  mode  control  scheme  has  been  developed 
for  the  9-node  plate  which  eliminates  spurious  modes  completely.  The  method 
is  based  on  retaining  the  formal  consistency  of  the  governing  equations  of  the 
systems,  which  is  equivalent  to  satisfying  the  patch  test.  As  a  result  of 
these  properties  of  the  spurious  mode  control  method,  results  which  have  been 
obtained  for  this  plate  show  nearly  the  optimal  h^  of  convergence. 

In  Chapter  1  of  this  report,  the  projection  methods  are  developed.  In 
the  first  section,  the  projection  method  is  developed  for  a  curved  beam  in 
order  to  illustrate  its  essential  features.  Section  2  then  uses  the  develop¬ 
ments  for  curved  beams  in  a  very  simple  fashion  to  develop  a  projection 
operator  for  the  3-node  triangular  plate  element;  this  element  uses  a  constant 
state  of  membrane  strain  and  either  the  discrete  Kirchhoff  theory  DKT  flexural 
element  or  a  C°  bending  formulation  with  one  quadrature  point.  Results  are 
then  given  in  Section  3  for  a  series  of  static  problems  and  two  dynamic 
problems  involving  the  collapse  of  shell  structures. 

In  Chapter  2,  the  procedure  for  spurious  zero  energy  mode  control  for  the 
9-node  element  is  developed.  The  development  is  first  given  in  the  setting  of 
the  Laplace  equation,  where  the  role  of  the  spurious  mode  control  procedure  is 
quite  transparent  and  the  important  role  of  maintaining  consistency  with  the 
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spurious  mode  control  procedure  is  illustrated.  Results  obtained  for  the 
Laplace  equation  show  that  by  using  the  spurious  mode  control  procedure,  con¬ 
vergence  in  the  L2  norm  of  order  h^  can  be  achieved  for  rectangular  meshes, 
whereas  some  deterioration  in  the  rate  of  convergence  results  when  the 
elements  are  curved  or  skewed.  However,  the  same  deterioration  in  the  rate  of 
convergence  occurs  if  selective  reduced  integration  is  used.  The  method  is 
then  applied  to  plate  problems.  A  large  series  of  plate  problems,  some  of 
them  involving  situations  where  convergence  is  very  slow,  such  as  a  rhombic 
plate  modeled  by  parallelogram  elements  are  considered.  In  all  cases,  very 
rapid  rates  of  convergence  are  achieved  with  this  element  and  no  spurious 
modes  are  present. 


STRAIN  PROJECTION  METHODS  FOR  CONSTANT  STRAIN  ELEMENTS 


1 .  Introduction 

The  development  of  an  effective  but  simple  shell  element  which  incorp¬ 
orates  the  effects  of  the  curvature  of  the  shell  but  avoids  locking  and 
spurious  kinematic  modes  is  essential  for  effective  and  economical  analysis  of 
shells  in  the  failure  domain.  However,  while  elements  such  as  the  16  node  C° 
element  are  quite  accurate,  their  complexity  and  high  cost  makes  them  unat¬ 
tractive  for  nonlinear  analysis.  Lower  order  or  flat  elements,  on  the  other 
hand,  tend  to  be  excessively  stiff  and  a  very  large  number  are  required  for 
accuracy. 

In  this  paper  we  sketch  the  development  of  a  simple,  curved,  triangular- 
shell  element.  The  element  has  the  following  advantages: 

1.  It  correctly  represents  rigid  body  motion. 

2.  It  correctly  represents  states  of  constant  membrane  strains  and 

constant  curvatures,  and  thus  allows  for  inextensional  bending  and 
eliminates  membrane  locking  [1]. 

3.  It  couples  bending  and  membrane  effects  within  an  element. 

4.  As  opposed  to  various  elements  based  on  selective  reduced  integ¬ 

ration,  it  possesses  no  kinematic  modes  [2-4]. 

5.  While  it  is  perhaps  the  simplest  curved  shell  element  (compare  [5- 
17])  the  element  yields  surprisingly  accurate  results,  often  superior 
to  those  obtained  by  more  complex  elements. 

The  basis  for  this  element  is  the  mode  decomposition  technique  described 
for  the  curved  beam  in  [18]  in  conjunction  with  a  shallow  shell  theory.  It 


should  be  noted  that  contrary  to  the  results  of  earlier  investigators  [12,19], 
use  of  a  shallow  shell  theory  in  our  case  does  not  diminish  the  element's 
performance  in  deep  shell  problems. 

The  mode  decomposition  technique  is  used  to  achieve  an  element  free  of 
membrane  locking  and  with  the  correct  relationship  between  membrane  and 
bending  effects.  The  DKT  (discrete  Kirchhoff  theory)  element  [20,21]  (see 
also  [22-24])  is  used  to  form  the  bending  part  of  the  stiffness  matrix.  This 
portion  of  the  stiffness  matrix  may  be  replaced  with  any  other  triangular 
plate-element  stiffness  matrix  provided  that  corners  are  the  only  nodal  points 
of  the  element.  However,  the  most  rigorous  justification  of  the  development 
presented  herein  is  related  to  the  DKT  element. 

To  make  the  paper  self-contained,  we  begin  with  a  short  presentation  of 
the  major  ideas  in  the  context  of  curved  beams.  This  is  followed  by  a 
development  of  the  triangular  shell  element.  Finally,  the  performance  of  this 
element  is  demonstrated  by  a  number  of  solutions  to  various  shell  problems. 
Some  general  remarks  conclude  the  paper. 

2.  Curved  Beam  Element 

A  conclusion  that  can  be  drawn  from  [25-32]  is  that  the  ability  to 
represent  independent  bending  and  membrane  strain  states  is  crucial  for  the 
success  of  a  curved  beam  or  shell  element.  However,  satisfaction  of  this 
requirement  is  difficult  and  has  only  been  accomplished  by  using  the  so-called 
assumed  strain  elements  [26,27,30].  Unfortunately,  this  approach  appears 
impossible  for  arbitrary  shell  elements  and  so  far  only  cylindrical  shell 
elements  have  been  formulated  by  this  method  [33-35].  Another  important 
conclusion  of  past  research  is  that  the  proper  inclusion  of  riyid  body  motion 


considerably  improves  element's  performance  [36-39],  but  this  is  not  central 
to  the  topic  of  this  paper. 

In  [18]  (see  also  [40])  a  different  approach  has  been  taken  to  meeting 
these  conditions.  It  allows  for  simultaneous  existence  of  both  membrane  and 
bendiny  strains  in  all  patterns  of  deformations  but,  at  the  same  time,  for  any 
given  set  of  nodal  degrees  of  freedom,  it  defines  certain  modes  of  deformation 
from  which  the  membrane  strain  energy  is  removed.  Since  only  the  bending 
strain  energy  is  assigned  to  these  modes,  we  call  them  bending  modes.  The 
remaining  portion  of  the  total  deformation  is  called  the  membrane  mode.  This 
modification  of  the  membrane  strain  energy  results  in  a  modified  and  better 
element  stiffness.  A  theoretical  justification  for  this  approach  to  curved 
beams  through  the  Hu-Washizu  variational  principle  is  given  in  [18], 

To  describe  how  the  bending  mode  and  consequently  the  membrane  mode  is 
defined,  consider  the  curved  beam  shown  in  Fig.  1.  To  properly  account  for 
rigid  body  motion  and  for  the  sake  of  simplicity,  we  consider  the  element  in  a 
corotating  frame  whose  x-axis  passes  through  both  of  its  ends.  A  theory  of 
shallow  structures  will  be  used  with  the  following  kinematical  relationships 

e  =  u,x  +  w°x  w,x  (la) 

<  =  -  w,xx  (lb) 

where  w°  describes  the  initial  shape  of  the  element,  u  and  w  are  displace¬ 
ments,  parallel  and  perpendicular  to  the  chord  respectively.  The  deform- 
ational  degrees  of  freedom  are 
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where  U21  is  the  relative  lonyitudal  displacement  of  the  right  end  with 
respect  to  the  left  one,  <t>j  are  the  deformational  rotations  and  4> j0"*'  the  total 
rotations  of  the  nodes. 

The  following  interpolations  are  used  for  u,  w  and  w° . 


u  ■  “21  ? 


w  =  L  +  (j^Ng) 

w°  =  L  (ajNj  +  <x2N2; 


where  a^,  are  rotations  associated  with  the  initial  shape  of  the  element,  L 
is  its  lenyth  and 
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In  view  of  Eq.  (lb)  it  is  evident  that  the  distribution  of  curvatures  is 
related  exclusively  to  w.  However,  based  on  Eqs.  (3a, 3b)  and  (la)  it  is  also 
evident  that  for  any  <  *  0,  when  wQ  *  0  then  e  *  0,  so  inextensible  bending 
states  are  impossible  for  the  approximation  given  in  Eq.  (3).  This  is 
precisely  the  cause  of  membrane  locking  discussed  in  [1].  However,  as  opposed 
to  the  reduced  integration  employed  in  [2],  membrane  locking  can  be  eliminated 
by  mode  decomposition  [18],  which  provides  a  more  rational  and  accurate 
method . 

In  the  decomposition  technique,  through  physical  arguments  certain  modes 
of  deformation  are  required  to  be  free  of  membrane  strain  energy.  An  operator 
is  then  developed  so  that  the  element  satisfies  these  requi rements .  For 
example,  in  a  curved  beam,  the  nodal  displacements  are  decomposed  into  bending 
and  membrane  modes.  In  bending  modes,  the  membrane  strain  energy  is  required 
to  vanish.  A  similar  approach  has  been  developed  to  avoid  shear  locking  [3]. 

To  define  the  bending  mode  of  this  decomposition,  note  that  in  inex- 
tensional  bending,  the  chord  of  a  curved  beam  must,  in  general,  change  in  the 
length.  This  longitudinal  displacement  u^  is  associated  with  the  rotations 
<t> ^  and  <t>2  and  define  the  bending  mode.  This  mode  of  deformation  must  be  free 
of  membrane  strain  energy  for  inextensional  bending  to  be  possible,  which  is 
imperative  if  membrane  locking  is  to  be  avoided. 

The  bending  mode  can  be  easily  determined  by  considering  that,  for  the 
displacement  fields  given  in  Eq.  (3b, 3c),  and  the  curvature  field  defined  in 
Eq.  (lb),  inextensional  bending  is  accompanied  by  the  following  extension  of 


the  chord 


(8) 


> 


and  E  is  Young  modulus,  I  the  moment  of  inertia  and  A  the  cross-sectional 
area.  The  element  stiffness  matrix,  as  usual,  is  obtained  from  Eq.  (8)  and 
given  in  terms  of  the  deformational  degrees  of  freedom  d  by 
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where 


C1  =  T5  “l  '  1U  a2 

r  _  2  1 

C2  TJ  2  U  1 


(10b) 

(10c) 

(10d) 


The  global  stiffness  matrix  is 


s  ■  (y2  -  y^/L 


c  =  (x2  -  x^/L 


(12b) 


It  is  important  to  note  that  inextensional  bending  is  now  possible  and 

A  A 

that,  through  the  nonzero  values  of  <i2  and  K13’  coupling  between  membrane  and 
flexural  behavior  is  introduced.  Hence  the  main  characteristic  of  a  curved 
element  is  retained. 

Remarks : 

1.  Although  a  physical  approach  has  been  presented  here,  a  formal 
equivalence  between  the  above  formulation  and  a  mixed  formulation  has 
been  established  in  [18];  the  technique  can  also  be  viewed  as  a 
strain  projection  method. 

2.  It  can  be  shown  that  this  formulation  is  also  equivalent  to  an 
"exact"  (without  any  modifications  or  reduced  integration) 
displacement  formulation  which  can  be  obtained  from  the  present  one 
by  replacing  the  linear  interpolation  for  u  given  in  Eq.  (3a)  by  a 
fifth-order  polynomial  that  satisfies  e  =  0  at  each  point  of  the 
beam. 

3.  The  same  approach  can  also  be  adopted  for  higher  order  elements, 

[18].  However,  the  higher  the  order  of  the  element,  the  less  severe 
the  membrane  locking  [2]. 


4.  For  a  straight  beam,  the  stiffness  matrix  obtained  from  Eq.  (8) 
reduces  to  two  matrices:  a  bar  matrix  and  a  beam  matrix  with  no 
interaction  between  them. 


Curved  Triangular  Shell  Element 


The  triangular  shell  element  has  three  nodes  at  its  corners  and  the 
corotational  plane  (x,  y)  is  always  defined  by  these  nodes  as  shown  in  Fig. 

2.  The  vectors  ej,  I  =  1  to  3,  are  unit  vectors  in  the  positive  directions  of 
the  local  coordinates  of  the  sides,  6  (U ,1 ) . 

The  nodal  degrees  of  freedom  at  each  node  are 


i\ot  -  (u  If.  u‘f,  w‘°*.  ,  .*")  I  -  1  to  3 


For  simplicity,  the  basic  developments  will  be  presented  in  terms  of 
deformation  nodal  displacements 


dT  =  (nT.  ST) 


n  =  (n1 ,  n2 .  n3) 


~  =  (9xl’  9yl  ’  9x2  *  9y2 ’  9x3’  9y3) 


(14a) 


(14b) 


(14c) 


where  nj  is  the  elongation  of  side  I  and  0xI,  6  j  are  the  deformation  nodal 
rotations  at  the  node  I,  which  are  given  by 
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8  .  =  0t°t  -  to 

xl  xl  x 


(15a) 


0  .  =  0*9*  -  to 

yi  yi  y 


(15b) 


where  tox  and  to^  are  the  rigid  body  rotations  of  the  element. 

The  membrane  state  of  strain  is  considered  constant  with  linear 
u  and  u  and  the  discrete  Kirchhoff  theory  (DKT)  is  used  to  describe  the 

*  J 

bending  properties  of  the  element.  In  the  DKT  element,  the  rotation  field  is 
quadratic  and  the  transverse  displacements  are  cubic  along  each  of  the  sides 
C20  j. 

It  should  be  noted  that  the  discrete  Kirchhoff  constraints  are  imposed  at 
three  points  along  each  side,  which  guarantees  orthogonality  of  the  tangential 
component  of  the  normals  to  the  midline  at  any  point  of  the  side  (see  [20,21] 
for  details).  Therefore  each  side  of  the  element  deforms  in  exactly  the  same 
fashion  as  the  curved  beam  based  on  Kirchhoff  (not  discrete  Kirchhoff)  theory, 
described  in  the  previous  section.  Thus,  if  the  in-plane  displacements  are 
linear  (strains  are  constant),  inextensional  bending  is  not  possible. 

Therefore,  a  modification  similar  to  that  employed  in  the  curved  beam 
element  is  made:  the  in-plane  elongations  which  should  accompany  inextensional 
bending  are  first  determined,  and  their  contributions  to  the  membrane 
strain  energy  then  are  removed. 

The  initial  curved  shape  of  the  element  is  described  by  w°(x,  y) ,  which 
along  each  of  the  sides  I  is  given  by 


V  V  *, 
W\ 


w°Uj)  =  a[  Ut)  +  a*  N2  Uj) 


.‘'v'vV'/v/vvvVvy''/'^  v 
v.va  .v  ;  -  .'a  . 


'•.’\-'vVVvVv\  s s'- y  - 
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where  aj  and  are  the  initial  rotations  of  the  midi i ne  for  side  I,  the 
length  of  side  I,  and  Nj(c)  are  given  by  Eqs.  (4). 

In  the  DKT  element,  the  transverse  displacements  along  each  side  are 
given  by 

w(£j)  =  tj  \  Ux)  +  <t>2  n2  )  (17 ) 


where 


>1  =  0xJeIy  "  9yJ  eIx 


(18a) 


♦2  =  9xl  eIy  '  9yl  elx 


(18b) 


and  0  =  I  -  1,  J  =  3  for  1-1=0. 

The  elongational  strain  and  curvature  along  each  edge  are  given  by 


e 


s 


u 


S.S 


+ 


(19a) 


-  w 


’ss 


(19b) 


From  the  similarity  of  Eqs.  (1)  and  (3)  to  (16),  (17)  and  (19),  it  is 
readily  apparent  that  each  side  of  the  plate  behaves  exactly  like  the  curved 
beam  described  in  the  previous  section.  Thus,  if  inextensional  bending  of  the 
shell  element  is  to  be  possible,  the  elongations  of  the  chords  associated  with 


bending  must  be  omitted  from  the  membrane  strain  energy. 

The  bending  elongation  component  of  each  side  is  given  by  the  counterpart 


of  Eq.  (5), 


nI  =  lU  *1  *1  (“4al  +  a2^  +  lu  *2  1 1  (“l  "  4a2^  (2°) 

Thus  if  the  total  strain  energy  is  taken  to  be 

u=7£TKbS+7(a-nb)  Km(a-nb)  (2i) 

then  bending  without  any  membrane  strain  energy  will  always  be  possible.  In 
the  above,  Kb  is  the  plate  bending  stiffness  of  the  element  expressed  in  terms 
of  its  6  deformation  modes,  e,  and  Km  the  membrane  stiffness  expressed  in 
terms  of  the  deformation  modes  n;  Km  is  identical  to  the  stiffness  of  the 
well-known  constant-strain  triangle  but  expressed  in  an  alternate  set  of 
degrees  of  freedom,  it  is  given  in  Appendix  A.  It  is  also  possible  to 
describe  K  directly  in  terms  of  nodal  displacements. 

It  is  convenient  to  express  nb  in  terms  of  e  by 

rjb  =  T  0  (22) 

where  T  is  given  in  Appendix  A.  Using  £q.  (21-22),  the  stiffness  of  the  shell 
element  can  be  shown  to  be  given  in  terms  of  the  deformation  degrees  of 
freedom  by 


K  = 


~m 

-TT  Km 


Km  T 
~m  ~ 


Kh  +  T  K„  T 
~b  ~  ~m  ~ 


(23) 


v 


In  the  above  equation,  the  effect  of  the  decomposition  becomes  transparent: 


1)  it  introduces  membrane-flexural  coupling  through  the  terms  K^T. 

2)  it  allows  the  membrane  stiffness  to  effect  the  bending  stiffness. 

These  effects  are  absent  in  any  flat  shell  element.  Thus,  the  most  important 

features  of  a  curved  element  are  embodied  in  this  element. 

Remarks: 

1.  In  this  element,  independent  states  of  arbitrary  curvatures  (bending 
modes)  and  membrane  strains  (membrane  modes)  are  possible.  Thus  the  ele¬ 
ment  can  always  bend  without  extension.  Although  not  all  shells  allow  for 
inextensional  bending  (see  [32]  for  example),  the  ability  of  the  element 
to  undergo  inextensional  bending,  regardless  of  its  shape,  is  beneficial, 
because  if  inextensional  bending  is  not  possible  in  a  particular  shell, 
the  appropriate  restrictions  on  the  bending  are  introduced  by  interactions 
between  elements,  but  not  the  element  itself. 

2.  For  linear  problems,  the  element  stiffness  matrix  of  Eq.  (23)  can  be 
expressed  in  closed  form,  so  no  numerical  integration  is  needed.  This 
results  from  the  fact  that  all  the  matrices  involved  in  Eq.  (23)  can  be 
expressed  in  a  closed  form  (see  [21]  and  [41]  for  details). 

3.  For  a  plate  problem  T  =  0  and  no  interaction  between  bending  and  membrane 
effects  exists. 

4.  Equation  (20)  is  fundamental  in  deriving  the  element  stiffness  matrix  for 
the  curved  triangular  shell  element.  It  is  valid  for  a  cubic  approxi¬ 
mation  for  the  initial  shape  of  the  element  sides  and  a  cubic  approxi¬ 
mation  of  transverse  displacements.  For  those  reasons,  the  DKT  element 
which  incorporates  cubic  transverse  displacements  along  each  of  its  sides 
seems  to  be  appropriate  for  the  present  development.  However,  other  ele- 


merits  can  be  formulated  with  cubic  interpolations  along  each  side.  The 
matrix  T  of  Eq.  (22)  would  then  remain  the  same  but  the  matrix  Kb  of  Eq. 
(23)  may  change.  Therefore  Kb  can  be  associated  with  other  formulations 
of  the  triangular  plate  element  without  changing  the  matrix  T.  In 
addition  to  the  results  obtained  with  Kb  associated  with  the  DKT  element, 
results  obtained  with  Kb  from  an  earlier  paper  [4]  will  be  presented  (see 
also  [3]).  This  matrix  was  found  to  closely  approximate  that  of  the  DKT 
element  although  it  is  simpler  and  more  economical  computationally. 


4.  Numerical  Results 

Results  are  presented  for  a  variety  of  linear  static  and  nonlinear  tran¬ 
sient  problems.  The  purpose  of  the  first  set  of  examples  is  to  illustrate  the 
enhanced  accuracy  of  the  curved  element  over  a  flat  element,  and  to  verify  the 
absence  of  membrane  locking.  Moreover,  some  deep  shell  problems  are  con¬ 
sidered  to  emphasize  that  elements  based  on  shallow  theory  do  converge  to  the 
deep  shell  solution.  The  second  set  of  problems  serves  simply  to  illustrate 
the  potential  of  this  element  in  highly  nonlinear  situations  such  as 
encountered  in  collapse  and  post-buckling  analysis. 

Results  are  presented  for  two  bending  stiffnesses:  the  DKT  stiffness 
[21]  and  the  C°  bending  stiffness  developed  in  [3,4], 

Bending  of  an  Infinite  Circular  Cylindrical  Shell 

The  problem  description  is  given  in  Fig.  3.  This  is  in  fact  an  arch  in  a 
state  of  plane  strain.  The  problem  has  been  solved  with  the  curved  shell 
element  and  with  a  flat  element.  The  results  normalized  with  respect  to  the 
analytic  solution  are  given  in  Table  1.  As  can  be  seen,  the  results  improve 


dramatically  when  the  initial  curvature  is  taken  into  account. 


For  a  different  arrangement  of  elements,  the  results  obtained  with  the 
flat  elements  may  be  better.  However,  in  the  analysis  of  general  shells,  the 
selection  of  a  mesh  is  not  as  easy  as  in  this  case.  We  therefore  have 
selected  the  mesh  shown  in  Fig.  3  to  illustrate  that  a  combination  of  flat 
elements  and  an  inappropriate  mesh  can  result  in  a  considerable  stiffening  of 
the  model.  Our  element  successfully  avoids  this  problem. 

Pinched  Cylinder  with  Free  Edges 

The  cylinder  and  the  related  data  are  given  in  Fig.  4.  This  problem  has 
been  analyzed  by  many  authors  and  their  results  are  compared  to  ours  in  Tables 
2  and  3.  The  total  number  of  degrees  of  freedom  was  the  same  for  all 
solutions.  Both  curved  and  flat  elements  with  both  the  DK.T  and  the  C°  bending 
formulations  were  used.  It  is  clear  that  our  element  compares  very  well  with 
the  most  accurate  elements  of  others,  including  formulations  which  were 
designed  specifically  for  cylindrical  shells,  [36],  [43]. 

It  is  interesting  to  see  that  in  this  case  the  improvement  achieved  with 
curved  elements  is  small.  This  is  perhaps  because  the  variation  of  displace¬ 
ments  along  the  meridian  of  the  shell  is  relatively  small  (in  spite  of  the 
concentrated  force)  due  to  the  small  ratio  of  length  to  radius  (see  Fig.  4). 
This  problem  is  not  a  demanding  test  on  the  performance  of  curved  shell 
elements. 

Pinched  Cylinder  with  Diaphragms 

The  problem  and  pertinent  data  is  depicted  on  Fig.  5.  The  results 
obtained  with  both  curved  and  flat  DKT  elements  along  with  those  presented  in 
[42]  are  reported  in  Table  4.  The  erratic  behavior  for  the  2x2  mesh  is  due 


to  the  inadequacy  of  the  shallow  shell  theory  for  the  relatively  deep  element 
emerging  from  such  a  coarse  mesh.  A  considerable  improvement  is  obtained  with 
curved  elements. 


This  problem  is  defined  in  Fig.  6  and  Table  5.  The  ends  of  the  panel  are 
simply-supported,  while  the  side  boundaries  are  fixed.  Both  material  and 
geometric  nonlinearities  were  considered.  The  load  was  applied  by  prescribing 
the  initial  velocity  given  in  Fig.  6  to  the  nodes  in  the  region  loaded  by  the 
explosive.  Figure  7  shows  an  undeformed  and  deformed  mesh.  Table  5  compares 
the  results  obtained  for  various  meshes  to  an  experimental  result.  It  can  be 
seen  that  the  convergence  to  the  experimental  value  is  relatively  slow.  The 
reason  for  this  however,  is  not  the  accuracy  of  the  element,  but  the  extreme 
localization  of  deformation  which  occurs  due  to  the  formation  of  plastic 
hinges. 


Collapse  of  a  Hollow  Column 


Figure  8  shows  the  simulation  of  a  hollow  column  loaded  axially.  The 
time  history  of  the  load  and  material  and  geometric  properties  may  be  found  in 
Fig.  3  and  Table  1  of  Ref.  [44],  respectively,  where  results  obtained  with  a  4 
node  quadrilateral  element  using  one-point  quadrature  and  hourglass  control 
are  also  given.  The  results  obtained  with  this  element  compare  well  with 
those  obtained  for  the  quadrilateral,  except  the  model  is  somewhat  stiffer. 
Note  the  severe  change  in  cross-section  which  accompanies  buckling. 
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APPENDIX  A 

The  matrices  Km  and  T  of  Eq.  (23)  will  be  presented  here  In  the  form 
consistent  with  the  deformational  degrees  of  freedom  given  in  Eqs.  (14b)  and 
(14c). 

To  obtain  the  required  form  of  Km  note  that  it  represents  membrane  strain 
energy  of  the  constant-strain,  flat  triangular  element 

U  ■  i  /  eT  D  £  dA  (24) 

a  ~ 

where  a  is  the  area  of  the  element 


D  = 


Eh 

(1-v2) 


1 

V 

0 


v  0 

1  0 


(25a) 


fex’  V  Yxy^ 


(25b) 


For  a  constant  strain  triangle  it  is  natural  to  replace  e  ,  e  ,  y  with 

a  y  xy 

strains  along  each  of  its  sides  Cj.  Since 

cT  -  -r-  ,  I  =  1.  2,  3  (26) 

l  *1 


and 


ex  eIx  +  £y  eIy  +  Yxy  eIx  eIy 


(27) 


we  obtain 


e  =  S'*1  ri 


where 


lelx 

*lely 

Zlelxely 

2 

2 

2e2x 

*2e2y 

*2e2xe2y 

2 

2 

3e3x 

*3e3y 

*3e3xe3y 

is  the  length  of  side  I  and  ejx,  e^  are  x  and  y  components  of  the  vector 
(Fiy.  2).  In  view  of  Eqs.  (24)  and  (28)  the  desired  form  of  the  matrix 


K  is 


*  A(S"1)  D  S'1 


To  obtain  the  matrix  T  of  Eq.  (22)  we  observe  that  it  defines  the 
elongation  nj  of  each  of  the  curved  sides  of  the  element  due  to  inextensional 


”*lC2ely  *'lC2elx 
-i2  Cle2y  *lCle2x 


'*2C2e2y 
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Deflection  Under  Point  load  for  Pinched  Cylinder 
with  End  Diaphragm* 


mesh 

on 

octant 

degrees 

of 

freedom 

OKT 

'  Dvorkin 

Curved 

flat 

3a the  C42] 

2  x  2 

22 

1 .324 

.  0.054 

■ 

4x4 

84 

0.777 

'  0.462 

- 

5  x  5 

130 

- 

1 

0.51 

6x6 

186 

0.269 

I  0.727 

1 

i 

8x3 

328 

0.948 

]  0.860 

i 

10  x  10 

• 

510 

0.988 

j  0.930 

j  0.83 

*  Results  are  normaltied  by  exact  solution,  0.12248  x  1C 


TABlE  5 

Material  and  geometric  constants  for  Impulsively 
loaded  cyl Indrical  panel 


Young's  modulus 
Density 

Poisson's  ratio 
Yield  stress 
Plastic  modulus 
Impulse  over  R. 


E  •  10.5  x  106 

p  »  2.5  x  ID*4  Ib-sec^/in* 

■j  •  0.31 

o  •  440C0.  psi 

Ep  •  o.  psi 

5650  In/sec 


TABLE  6 


Final  Displacements  of  the  Cylindrical  Panel 
for  Various  Meshes  »ith  Ilyushin  field  Condition 


Mesh  for 

Ha)  f*pdnel 

Displacement 
at  y  •  5.2S 

Displacement 
at  y  *  9.4’ 

CP'J  tim 

I  DM  3033 

6  x  16 

0.317 

0.40) 

33 

3  i  16 

1  .043 

0.448 

138 

10  x  20 

1 .081 

0.462 

162 

12  «  24 

1 .124 

0.473 

291 

16  x  32  5 

1.183 

0.530 

570 

experimental  [ij] 

1.28 

0.70 

Curved  bean  nonenclature 


R'4.833  lh*0.094 

L- 10.35  „  |P- 100.0 

E*10.5*10°  1  jt)*0. 01548 

0*0.3125  P  lP-0.1 


Figure  1.  Example  2.  Pinched  cylinder  with  free  edges. 
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Figure  5.  Example  3.  Pinched  cylinder  with  end  diaphragms. 


CHAPTER  2 

SINGULAR  MODE  CONTROL  IN  9  NODE  PLATE  ELMENT 


1.  Introduction 

The  numerical  quadrature  of  higher  order  elements  is  a  critical  issue 
which  must  be  treated  with  considerable  care.  It  was  recognized  by 
Zienkiewicz  et  al  [1]  that  in  order  to  obtain  reasonable  rates  of  convergence 
for  C°  plate  and  shell  elements,  reduced  integration  of  the  shear  terms  is  a 
necessity.  Recently  it  has  become  clear  that  for  curved  elements, 
overintegration  of  the  membrane  terms  [2-5]  also  impairs  convergence;  this 
phenomenon  has  been  termed  membrane  locking.  In  fact,  this  phenomenon  also 
occurs  in  the  application  of  continuum  elements  to  curved  members  [5], 

On  the  other  hand,  reduced  integration  leads  to  the  appearance  of 
spurious  singular  modes,  and  several  of  the  elements  which  have  appeared  in 
the  literature  for  plate  analysis  cannot  treat  the  modal  analysis  of  unsup¬ 
ported  plates  because  they  rely  on  boundary  conditions  to  suppress  these 
spurious  modes;  see  for  example  the  results  in  [6]  on  the  plate  element  of 
[7j.  These  modes  are  also  known  as  kinematic  modes,  spurious  zero-energy 
modes,  and  for  the  4-node  quadrilateral,  as  hourglass  or  keystone  modes. 

In  the  4  node  continuum  and  plate  elements,  considerable  progress  has 
recently  been  achieved  in  controlling  hourglass  modes.  Kosloff  and  Frazier 
[8]  first  recognized  that  in  controlling  spurious  singular  modes,  care  must  be 
taken  to  ensure  that  the  resulting  element  will  pass  the  patch  test  and  not 
adversely  affect  rigid  body  and  constant  strain  patterns.  While  the  method  of 
[8]  involved  the  solution  of  systems  of  equations,  in  [9-10]  explicit  forms  of 
the  operators  for  controlling  singular  modes  were  developed.  In  [11],  these 
forms  were  derived  from  the  consistency  conditions  (which  are  equivalent  to 
the  patch  test  [12]).  This,  in  constrast  to  [8],  makes  the  development  of 
simple,  explicit  forms  of  the  spurious  mode  control  feasible  for  a  large  class 


of  elements;  by  usiny  a  mixed  variational  principle,  the  identification  of 
control  parameters  for  nonlinear  materials  was  also  achieved.  Applications  of 
these  concepts  to  the  diffusion  equation,  plate,  and  shell  problems  are  yiven 
in  [12-14],  respectively.  The  development  of  plate  elements  of  proper  rank 
has  also  been  undertaken  via  hybrid  [15]  and  mixed  formulations  [16], 

In  this  paper,  the  consistency  conditions  will  be  used  to  derive  a 
spurious  mode  control  for  the  9  node  Lagranye  element.  Pawsey  [17]  evidently 
first  noted  their  appearance  in  higher  order  elements.  Shortly  thereafter. 
Cook  [18]  proposed  a  control  which  consisted  of  simply  addiny  a  spring  between 
the  center  nodes  of  all  elements  and  a  fixed  point.  This  method  obviously 
would  not  pass  the  patch  test,  but  its  performance  in  simple  problems  is  not 
bad.  The  presence  of  spurious  singular  modes  in  the  9-node  Lagrange  elements 
with  2x2  quadrature  was  demonstrated  numerically  in  continuum  elements  and 
the  diffusion  equation  in  [19]  and  [20],  respectively.  Cook  and  Zhao-hua  [21] 
have  proposed  a  control  method  which  consists  of  perturbing  the  9-node  stiff¬ 
ness  by  that  of  the  8-node  serendipity  element.  This  element  passes  the  patch 
test  but  requires  condensation  of  the  interior  node,  which  is  awkward  for 
nonlinear  applications. 

In  this  paper,  a  method  for  controlling  (or  stabilizing)  the  spurious 
singular  modes  by  using  the  consistency  conditions  in  the  manner  of  [11]  is 
developed.  The  control  is  first  developed  for  the  diffusion  operator  because 
it  is  the  simplest  setting  in  which  these  modes  appear;  Section  2  gives  the 
discrete  form  of  the  diffusion  equation  to  be  used  here.  Section  3  describes 
the  spurious  mode  control  procedure  and  Section  4  gives  results  for  selected 
problems.  The  important  feature  of  these  results  is  that  with  2x2  quad¬ 
rature  of  the  9-node  element  and  hourglass  control,  the  convergence  rate  of 


where  upper  case  subscripts  refer  to  node  numbers  and  n^  is  the  number  of 
nodes  possessed  by  the  element.  Followiny  standard  finite  element  procedures, 
the  conductance  matrix  can  be  shown  to  be  given  by 


r  -  -  l  "i.i  *u  "j.j  d(! 


(2.3) 


The  above  is  often  written  in  the  form 


Ke  =  /  BT  D  B  dfl 


(2.4) 


where 


B  = 


=  (NI’X 

u 

1  NI.,* 

(2.5a) 


D  = 


“ll  “12 

“12  “22 


(2.5b) 


where  the  symmetry  of  the  coefficients  has  been  invoked;  B  is  called  the 

discrete  form  of  the  gradient  operator  since  (see  Eqs.  (2.2)  and  (2.5a)) 


n. 


'N 

[ 

1  =  1 


u’i  =  l  6il  UI  ~i  ~ 


e  =  bj  ue 


(2.6) 


and  the  transpose  of  B  is  the  discrete  form  of  the  divergence  operator. 
For  the  9-node  Lagrange  element,  the  shape  functions  Nj  are 


Nj  =  JCj(C)  *Kln)  I  =  3 ( J -1 )  +  K 


(2.7) 


where  jU)  are  the  standard  one-dimensional  Lagrange  interpolation  functions 
in  terms  of  the  reference  coordinates  (see  Fig.  1)  over  the  points  (-1,  0, 
+1)  which  are 


siCU-l) 

2 


c£2U)  =  1  -  K2 


(2.8) 


=  -  5  (5+D 
J  2 


The  reference  plane  is  related  to  the  physical  plane  by 


xi  3  j l=1  NI  xil  3  ^*i 


(2.9) 


where  are  the  coordinates  of  node  I;  the  nodal  coordinates  are  sometimes 
arranged  in  2  column  matrices  x^ ,  where  x^  =  x,  x2  =  y. 

When  numerical  quadrature  is  used,  the  conductance  matrix,  Eq.  (4),  is 
evaluated  by 


Ke  =  l  p  BT  (xQ)  D  (x  )  B  (xQ) 
a=I 


(2.10) 


where  x^  =  (x^,  y^),  a  =  1  to  nQ,  are  the  quadrature  (integration)  points 


and  the  quadrature  weights  including  the  Jacobian.  The  above  can  be 
written  as 


I  a( 
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The  rank  of  Ke  is  at  most  equal  to  the  rank  of  B.  This  is  a  consequence 
of  the  fact  that  the  dimension  of  the  null-space  of  §  is  equal  to  that  of  the 
null-space  of  Ke  because  if  5  u  =  0,  then  Ke  u  =  0;  see  Eq.  (2.11).  With  this 
quadrature  scheme  the  rank  of  §  could  be  at  most  B,  since  B  has  8  rows; 
however  one  of  these  8  rows  is  a  linear  combination  of  the  other  7,  so  the 
rank  of  B  and  Ke  is  7;  this  rank  deficiency  has  also  been  reported  in  [20]. 


3.  Spurious  Mode  Control  for  Laplace  Equation 

The  essential  feature  of  an  effective  control  of  spurious  singular  modes 
is  that  the  stabilization  matrix  should  not  result  in  a  violation  of  the 
consistency  conditions.  The  consistency  conditions,  which  are  in  a  sense 
equivalent  to  the  patch  test  [22,  23]  require  that 

£i0*j=«ij  for  all  a  (3.1) 


bls=0  or  S  s  =  0 


(3.2) 


S  =  [1,  1,  1,  1,  1,  1,  1,  1,  1] 


(3.3) 


Equations  (3.1)  and  (3.2)  were  verified  algebraically  for  rectanyular  elements 
and  numerically  for  a  large  variety  of  curved  elements. 

Note  that  Eq.  (3.1)  corresponds  to  the  requirement  that  (cf. 

Eqs.  (2.5a),  (2.9)),  i .e.  that  the  gradient  of  a  linear  field  be  computed 
exactly,  while  Eq.  (3.2)  implies  that  any  constant  field  has  a  zero 
gradient.  Using  Eq.  (2.11),  Eq.  (2.16)  can  be  seen  to  imply  that 


Ke  s  =  0 


(3.4) 


Thus  s  is  in  the  null -space  of  B  and  Ke.  The  vector  space  spanned  by  s  is 
called  the  proper  null-space  of  Ke  because  the  gradient  is  expected  to  vanish 


for  a  constant  field  u(x,y)  which  is  associated  with  the  nodal  values  ue  =  s. 


It  can  also  be  shown  that 


b.  h  =  0  or  B  h  =  0 


(3.5) 


K  h  =  0 


(3.6) 


where 


h  =  [+1,  -1,  +1,  -1,  0,  -1,  +1,  -1,  +1] 


(3.7) 


Equation  (3.5)  has  been  verified  algebraically  for  rectangular  elements  and 


checked  numerically  for  a  large  number  of  rhombic  and  curved  elements;  Eq. 
(3.6)  follows  from  (3.5)  because  of  Eq.  (2.11).  The  one-dimensional  vector 
space  of  h  is  called  the  improper  null -space  of  B  because  the  gradient  of  the 
function  u  associated  with  the  nodal  values  ue  =  h  is  obviously  not  zero,  as 
can  be  seen  from  a  depiction  of  this  mode  in  Fig.  2. 

Simple  examination  of  Eqs.  (3.3)  and  (3.7)  reveals  that  the  vector  h  is 
also  orthogonal  to  s,  i.e. 


hTs  =  0 


(3.8) 


so  that  the  entire  complement  to  the  proper  null  spaa.  of  s  is  spanned  by  the 
9  vectors  b.  and  h. 

The  control  of  the  improper  (spurious)  mode  is  accomplished  by  defining 
an  additional  generalized  gradient  g  by 


9  =  iV 


and  a  generalized  flux  by 


q  =  -  a  g 


(3.9a) 


(3.9b) 


where  the  determination  of  a  is  described  later. 

The  resulting  element  conductance  matrix  is  then  given  by 


Ke=BTDB+ayyT 


(3.10) 


\  • 
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where  the  second  term  in  the  above  is  the  stabilization  matrix  [6,12]. 

If  the  above  stiffness  matrix  is  to  meet  the  consistency  requirements  on 
linear  fields,  it  is  sufficient  and  necessary  that  g  vanish  for  any  nodal 
values  associated  with  the  linear  field  [11] 

u  =  c0  +  ciX  +  c2y  =  cQ  +  c^  ■  (3.11a) 

that  is,  for  nodal  values  of  u  given  by 

u  =  cQ  s  +  ci  xi  (3.11b) 

the  generalized  gradient  g  must  vanish.  This  is  equivalent  to  the  requirement 
that 

IT  (c0  2  +  ci  *,)  =  0  t3-1*) 

for  all  c-j ,  i  =  0  to  2.  Since  the  10  vectors  bia,  s  and  h  span  the  9 
dimensional  vector  space  R9,  the  most  general  form  of  i  is  given  by 

4 

*  =  l.  [ala  fela+  32a  ha*  +  ^  +  a10§  (3‘13) 

a=l 

Substituting  the  above  into  Eq.  (3.12)  and  using  Eqs.  (3.1)  (3.2),  (3.5)  and 
(3.8)  yields 


(3.14) 


co(alU  ~  ^  +  cl^a9~  ~  +  alU^  i  +  l  aia  ) 


o*l 


c2(a9^T^  +  al0§^  +  I  a2a}  =  0 

a=l 


Since  the  above  must  hold  for  all  c-j,  it  follows  that 


(3.15a) 


a9  +  l.  alc  =  0 


o=l 

4 


a9  JTi!  *  I  >2a  •  ° 

0=1 


(3.15b) 

(3.15c) 


Equations  (3.15)  represent  2  conditions  on  the  8  parameters  aia.  Therefore 
considerable  leeway  is  available  in  the  choice  of  these  constants  if 
consistency  is  only  required  with  respect  to  linear  fields.  It  can  be  met  by 
letting 


ai.  ’  -  7  asiTi- 


a  =  1  to  4 


a2«  ■  -  7  a9SV 


a  =  1  to  4 


which  with  Eqs.  (3.13)  and  (3.15)  gives 


1  =  £  ‘  -  C(i3Tii)  l  Sia  +  l*2 J 


_.~2a 


(3.16a) 

(3.16b) 


(3.17) 


The  common  constant  ag  has  been  omitted. 

For  purposes  of  constructing  the  stabilization  matrix,  (3.10),  the 
constant  a  is  expressed  in  terms  of  a  stabilization  parameter 
e  by 


ea 


a  = 


4 

I 

200  a=l 


11 


P  b  b. 
a  ~la  ~1a 


(3.18) 


In  this  form  e  =  1  closely  approximates  the  exactly  integrated  matrix  for 
rectangles  when  a  is  isotropic. 

Remark .  This  procedure  involves  8  arbitrary  constants,  a^a  and  a2a,  so  it 
would  be  advantageous  if  consistency  of  the  quadratic  and  biquadratic  poly¬ 
nomials  could  also  be  enforced.  However,  the  resulting  consistency  equations 
are  singular  and  we  have  not  been  able  to  circumvent  this  difficulty. 
Furthermore,  quadratic  consistency  is  not  satisfied  by  the  9-node  element  when 
the  element  is  not  a  parallellogram,  that  is, 


(3.19) 


when  u  is  a  quadratic  function  of  (x,y)  and  u  are  its  nodal  values. 

Remark:  Equation  (3.17)  can  also  be  written  in  a  form  easily  implemented  on  a 
computer  as  follows 


(3.20) 


and  x’  are  the  4  quadrature  points. 


4.  Numerical  Examples  for  the  Laplace  Equation 

Several  problems  were  solved  for  the  purpose  of  examining  the  performance 
of  this  spurious  mode  control  procedure,  with  particular  emphasis  on  the  rates 
of  convergence  and  the  effects  of  varying  the  parameter  e.  All  of  the 
problems  solved  were  linear  and,  except  for  one,  steady-state.  For 
simplicity,  the  dependent  variable  u  is  called  the  temperature  throughout  this 
section. 

Example  1  is  defined  in  Fig.  3.  Three  meshes  were  used  for  this 
problem.  The  temperature  distribution  along  the  axis  B-C  is  compared  with  the 
analytic  solution  in  Fig.  4.  Figure  5  shows  the  results  for  the  fine  mesh 
with  e  =  0.1  and  1.0,  respectively.  It  can  be  seen  that  in  this  problem  the 
effect  of  the  stabilization  parameter  e  is  minimal.  Because  this  problem  has 
a  prescribed  temperature  along  all  boundaries  (Dirichlet  type),  2x2 
quadrature  without  stabilization  does  not  lead  to  any  spurious  modes,  so  these 
solutions  primarily  show  that  the  introduction  of  this  spurious  mode  control 
procedure  does  not  distort  the  solution. 

Example  2  is  described  in  Fig.  6.  Again,  three  meshes  were  used  as  shown 
in  Fig.  6.  Four-fold  symmetry  was  used  in  these  problems  and  the  lines  of 
symmetry  are  considered  insulated,  i.e.  the  gradient  vanishes  along  the  lines 
x  =  0  and  y  =  0.  Since  this  is  a  natural  boundary  condition,  it  need  not  be 
explicitly  enforced.  The  solutions  for  three  different  values  of  e  are  given 


in  Fig.  7.  This  problem,  because  of  the  shortage  of  essential  boundary 
conditions,  can  exhibit  spurious  modes,  and  as  can  be  seen,  for  the  lowest 
values  of  e,  oscillations  in  the  temperature  u  are  quite  clear.  For 
e  =  l.U,  the  oscillations  are  almost  completely  eliminated  and  the  solution 
agrees  well  with  the  analytic  solution.  Note  that  the  results  with 
insufficient  stabilization  do  not  oscillate  about  the  correct  solution. 

Example  3  is  identical  to  example  2  except  that  a  Dirichlet  boundary 
condition  is  applied  to  the  circumference,  see  Fig.  6.  Results  for  the  three 
meshes  in  Fig.  6  are  given  in  Fig.  3,  and  compared  to  the  analytic  solution; 
results  for  two  values  of  e  are  given  in  Fig.  9.  Both  agree  well  with  the 
analytic  solution. 

Figure  10  shows  the  1_2  norm  of  error  as  a  function  of  element  size  h  for 
examples  1  and  3,  which  are  rectangular  and  circular  domains,  respectively. 
The  L9  norm  was  computed  by  evaluating 


E2  -  f  f  (uFEM  -  uandl*tic)  « 
e=l 


and  5x5  quadrature  was  used  to  evaluate  the  above  integral  in  each  element. 
Results  for  both  the  stabilization  procedure  with  e  =  1  and  3  x  3  quadrature 
are  given.  For  the  square  domain,  where  the  elements  are  all  straight  sided, 
the  slope  of  this  curve  is  2.95  for  both  the  stabilization  procedure  and  3x3 
quadrature.  Thus,  for  this  case  the  stabilization  procedure  has  no  undesir¬ 
able  effects  on  the  rate  of  convergence.  In  the  circular  domain,  example  3, 
where  some  of  the  elements  are  curved  and  irregular,  the  situation  is  more 
complex.  For  the  meshes  studied,  the  absolute  accuracy  of  3  x  3  quadrature  is 
better  than  that  of  stabilization  procedure.  However,  3x3  quadrature  exhi- 
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bits  a  definite  reduction  in  the  rate  of  convergence  with  decreasing  element 
size  which  is  absent  in  the  stabilized  element.  Average  slopes  are  1.9  for 
2x2  with  stabilization,  1.4  for  3x3  quadrature,  so  for  these  irregular 
meshes,  neither  3x3  quadrature  nor  the  stabilization  procedure  retains  the 

3 

optimal  h  convergence  of  this  element.  The  results  for  even  higher  order 
quadrature,  4x4  and  5x5,  are  almost  identical  to  that  for  3x3 
quadrature. 

Figure  11  is  included  to  show  the  effect  of  spurious  modes  more 
dramatically.  In  this  problem,  a  circular  domain  with  an  insulated  outer 
boundary  was  considered.  A  point  source  which  is  a  step  function  in  time  is 
applied  at  the  center  of  the  domain.  The  time  step  in  this  problem  is  one, 
and  evidence  of  spurious  oscillations  is  quite  clear  within  20  time  steps. 
Within  40  time  steps,  the  entire  solution  is  dominated  by  the  spurious 
oscillations.  Note  that  in  this  plot  all  nodes  are  connected  by  solid  lines, 
so  each  9-node  element  is  represented  by  4  adjacent  elements  in  Fig.  11. 

5.  Discrete  C°  Plate  Equations 

The  general  theory  for  C°  plates,  which  are  often  called  Mindlin  plates, 
is  presented  in  [24],  In  this  theory,  the  deformation  of  the  plate  is 
described  by  three  dependent  variables;  a  transverse  deflection  w(x,y)  and 
rotations  e  (x,y)  and  8  (x,y).  The  element  has  9  nodes  with  3  degrees  of 

a  jr 

freedom  per  node  and  the  three  fields  are  approximated  in  the  element  by  the  9 
shape  functions,  Eqs.  (2.7  -  2.8)  in  the  form 

nN 

w(x,y)  =  l  Ntwt  =  N  w  (5.1) 
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The  element  stiffness  matrix  Ke  is  given  by 


5e-  /Ae  sj  Sb  8b  *  /AeSTsSsSsdA 


(5.2a) 

(5.2b) 


(5.3) 


bending  stiffness 


shear  stiffness 


where  Ae  is  the  area  of  element  e,  and 
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Here  h  is  the  thickness,  E  is  Young  s  modulus,  G  is  the  shear  modulus, 
k  is  the  shear  factor  and  v  is  Poisson's  ratio. 

This  element  locks  if  3  x  3  quadrature  is  used  for  both  the  bending  and 
shear  stiffnesses;  therefore  selective  reduced  integration  is  recommended 
[24].  In  this  paper,  we  will  use  2x2  quadrature  for  both  the  shear  and 
bending  stiffness  and  stabilize  the  kinematic  modes. 

For  the  purpose  of  analyzing  this  stiffness,  we  perform  the  steps 
indicated  in  Section  2  in  going  from  Eq.  (2.4)  to  Eq.  (2.10)  and  write  the 
element  stiffness  for  2x2  quadrature  in  the  form 


Ke  =  B7  D  B 


(5.8) 


where 
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where  xjjj,  a  =  1  to  n^,  are  the  Gauss  quadrature  points;  in  this  case  ng  =  4. 

The  number  of  rows  in  the  matrix  B  is  20  so  the  rank  of  Ke  is  at  most  20; 
the  number  of  degrees  of  freedom  is  27;  (9  nodes  x  3  degrees  of  freedom/per 
node).  Thus,  since  there  are  3  rigid  body  modes,  the  rank  deficiency  of  the 
element  stiffness  is  at  least  4.  The  singular  modes  are  given  in  Table  1. 


6.  Spurious  Mode  Control  for  Plate  Element 


The  spurious  modes  which  are  of  principal  concern  are  the  w,  9  and  e 

X  y 

modes.  The  process  of  control  is  analogous  to  that  in  the  Laplace  equation, 
except  that  whereas  the  singular  modes  in  the  Laplace  equation  can  be 
controlled  quite  stiffly  with  values  of  t  of  the  order  of  1,  some  of  the  modes 
in  the  plate  problem  are  associated  with  locking  and  therefore  care  must  be 
taken  in  their  control. 

For  the  purpose  of  the  control  of  these  3  modes,  three  generalized 
strains  are  defined  by 


The  structure  of  the  proper  null -space  for  the  plate  is  more  complex  than 


I 


t: 


for  the  diffusion  equation;  it  is  spanned  by  the  3  vectors  which  are  called 
riyid  body  modes  in  Table  1.  These  vectors  include  both  the  linear  field 
vectors  and  the  constant  vector  s.  In  spite  of  this  added  complexity,  the 
basic  conditions  on  x<  can  be  reduced  to  those  on  x  in  the  diffusion  equation: 

1.  ^  should  not  affect  linear  fields: 

2.  should  span  the  complement  of  the  proper  null-space. 

Therefore,  the  same  procedure  as  given  in  Section  3  can  be  used  to 

develop  and  it  follows  that  all  three  of  the  x^  operators  are  identical  and 
given  by  Eq.  (3.17) . 

The  generalized  stresses  are  given  by 


Q.  =  c"  q. 

1  Hl 


for  i  =  1  to  3, 
no  sum  on  i 


(6.2) 


where  c?  are  given  by 


„H  rwh  HD 
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(6.3a) 


c2  =  C3  =  reH  <Gh 


(6.3b) 
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The  constant  H  was  herein  obtained  by  2  x  2  quadrature,  but  it  is  likely  that 
estimates  of  H  with  sufficient  accuracy  could  more  easily  be  obtained. 

The  generalized  strain  qw,  as  in  the  4  node  quadrilateral  [13],  is 
associated  with  locking,  so  rw  should  be  small;  suggested  values  are: 

0.01  <  r  <  0.1  .  In  this  range,  the  results  have  been  found  to  be  almost 
independent  of  rw.  The  normalization  of  Eqs.  (6.3)  allows  rw  to  be  applicable 
to  elements  with  a  wide  variety  of  shapes  and  aspect  ratios,  for  we  have  found 
no  locking  or  evidence  of  spurious  inodes  with  rw  in  this  range.  The 
generalized  strains  q2  and  qj  are  not  associated  with  a  locking  mode,  so  r0  is 
usually  chosen  of  order  1.0. 

The  element  stiffness  matrix  is  given  by 


K  =  K<2x2>  +  KH 


(6.5a) 


KH  = 


c?XIT 


0 

0 


0 

H  T 
c2  X  X 

0 


0 

0 

_H  T 
c3  X  X 


(6.5b) 


H  (2x2)  * 

where  K  is  the  stabilization  matrix  and  K'  the  standard  element  matrix 

t+m  /"W 

obtained  by  2x2  quadrature. 

If  the  degrees  of  freedom  of  the  element  are  arranged  in  the  conventional 
order  with  all  degrees  of  freedom  at  each  node  in  sequence,  i.e.  with 
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then  K  can  be  easily  computed  by  the  following  formula 

For  1=1  to  9, J=1  to  9 

C1  YI  YJ  L  =  31  -  2,  K  =  3J  -  2 

C2  YI  YJ  L  =  31  -  1,  K  =  3J  -  1  (6.7) 

C3  YI  YJ  L  -  31,  K  =  3J 

U 

All  other  are  zero. 

Thus  the  implementation  of  this  element  involves  simply  the  standard 
computation  of  the  element  stiffness  by  2  x  2  quadrature,  followed  by  the 

Li 

computation  of  x  by  Eq.  (3.17)  (or  (3.20))  and  K  by  the  above. 

7.  Numerical  Examples  for  Plate 

Several  problems  were  solved  to  examine  the  performance  of  this  spurious 
mode  control  procedure  in  the  plate  element. 

The  first  problem,  example  5,  is  a  clamped,  circular  plate  of  radius  R 
subjected  to  a  uniform  load  q.  Results  obtained  Mesh  B  from  Fig.  6  for  two 
values  of  rw  are  compared  to  the  analytic  solution  in  Fig.  12.  In  both  cases, 
the  e  -  hourglass  control,  rQ  =  1.0.  The  results,  as  can  be  seen,  are 
independent  of  r  in  this  range,  and  agree  very  well  with  the  analytic 
solution. 

Figure  13  shows  the  results  for  the  same  problem  with  rw  =  0.1  for  the  3 
meshes  shown  in  Fig.  6.  Even  with  the  coarse  mesh,  the  results  agree  well 


Figure  14  shows  the  convergence  rates  for  this  circular  clamped  plate 
problem  and  a  uniformly- loaded  square  plate  problem  with  clamped  supports; 
the  meshes  shown  in  Fig.  3  were  used  for  a  quarter  of  the  plate  in  the  square 
plate  problem.  The  L2  error  is  here  defined  as 

E2  -  /  (w™  -  wana'»t1c)2«  (7.1) 

e=i  ae 

and  5x5  quadrature  was  used  to  evaluate  the  above  integral  in  each 
element.  In  all  cases,  rw  =  0.1,  rQ  =  1.0  . 

Several  points  are  of  interest: 

i)  The  rate  of  convergence  in  w  for  the  stabilization  method  for  the  square 
plate  is  almost  3.0. 

ii)  The  rate  of  convergence  for  the  circular  plate  problem  with  stabilization 
is  2.6.  While  the  initial  rate  with  select! ve-reduced  integration  is  also 
2.6,  the  convergence  rate  diminishes  as  the  mesh  is  refined. 

iii)  The  behavior  with  an  inconsistent  hourglass  control  x  =  J2»  i»e. 
where  the  last  term  in  Eq.  (3.17)  is  omitted,  is  similar  to  that 
with  selective- reduced  integration. 

While  the  improved  rate  of  convergence  for  the  circular  plate  as  compared 
to  the  Laplacian  on  a  circular  domain  may  be  puzzling  at  first,  it  is  probably 
attributable  to  the  omission  of  the  rotations  from  the  error  in  Eq.  (7.1). 

The  performance  of  the  method  in  a  problem  characterized  by  severe 
singular  modes  is  shown  in  Figs.  15  and  16,  which  is  a  uniformly-loaded  square 
plate  with  corner  supports,  such  as  that  studied  in  [13],  For  small  values 
of  rw,  spatial  oscillations  in  w  are  clearly  evident.  Figure  16  shows  the 
normalized  displacement  of  a  node  next  to  the  center  for  the  corner  supported 


and  clamped  square  plates.  Note  that  as  for  the  4-node  quadrilateral  [13], 

the  stabilization  procedure  gives  acceptable  results  for  a  large  range 

of  r  ,  so  the  results  are  not  sensitive  to  its  selection.  The  clamped  plate 

locks  for  large  values  of  r  but  exhibits  no  modes  as  r  tends  to  zero, 

w  w 

whereas  the  corner  supported  plate  diverges  for  small  r  but  does  not  lock. 

w 

Only  the  intermediate  values  give  acceptable  solutions  to  both  problems. 

Figure  17  shows  the  displacements  for  a  circular  thick-plate  subjected  to 
a  point  load  at  its  center.  The  boundaries  are  clamped  and  the  following 
parameters  were  considered:  E  =  1.09  x  106  psi;  v  =  0.3;  thickness  h  =  2  in; 
radius  R  =  5  in.  Mesh  8  in  Fig.  6  was  used. 

In  this  case,  2x2  quadrature  leads  to  near  singularity  of  the  assembled 
stiffness,  and  the  results  with  no  stabilization  exhibit  marked  oscil¬ 
lations.  Effective  suppression  of  oscillations  in  this  case  requires  a  larger 
value  of  rw  (0.1)  than  in  any  other  problem  we  have  solved,  and  the  displace¬ 
ment  at  the  center  is  more  sensitive  to  r  . 

w 

Figure  18  shows  the  results  obtained  for  this  element  in  the  well-known 
"single-element  twist"  problem  [25]. 

The  last  well-known  difficult  problem  is  the  rhombic  plate.  The  two 
meshes  shown  in  Fig.  19  were  considered.  Displacements  and  moments  are 


reported  in  Figs.  20  and  21. 
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1.  Nomeclature,  local  node  numbers  and  quadrature  points  for  9-node  Lagrange 
element. 

2.  Spurious  singular  mode  h. 

3.  Problem  description  for  example  1:  square  domain  with  prescribed 
temperature  along  boundary. 

4.  Temperature  u  along  B-C  for  the  three  meshes  for  example  1  with  e  =  1.0. 

5.  Temperature  u  along  B-C  for  example  1  with  medium  mesh 
and  e  -  0.1  and  1.0. 

6.  Problem  description  for  example  2  and  3:  circular  domain  with  Neumann  and 
Dirichlet  boundaries. 

7.  Temperature  u  along  the  radius  coincident  with  the  x-axis  for  the  medium 
mesh  e  =  10“  ,  10”  ,  and  1.0  for  example  2. 

8.  Temperature  u  along  the  radius  coincident  with  the  x-axis  for  the  three 
meshes  for  example  3  with  e  =  1.0. 

9.  Temperature  u  along  the  radius  coincident  with  the  x-axis  for  example  3 
with  e  =  0.1  and  1.0. 

10.  Convergence  rates  for  examples  1  (square)  and  3  (circular)  with  2x2 
quadrature,  e  =  1.0  and  3  x  3  quadrature. 

11.  Temperature  distributions  for  example  4,  showing  the  evolution  of  spurious 
oscillations  with  2x2  quadrature  in  the  absence  of  stabilization. 

12.  Displacement  w  along  a  radius,  line  B-C  in  Fig.  6,  for  a  clamped,  circular 

plate  for  two  values  of  r,,. 
r  w 

13.  Displacement  w  along  a  radius,  line  B-C  in  Fig.  6,  for  three  meshes  for 
the  clamped,  circular  plate. 

14.  Convergence  in  L?  norm,  Eq.  (7.1),  for  clamped  square  and  circular  plates 
with  stabilization  (r  =  0.1,  rQ  =  1.0  )  and  for  selective  reduced 
integration;  the  convergence  rate  when  ^  =  h  is  also  shown. 

15.  Displacements  along  a  centerline  for  a  uniformly  loaded  a  x  a  square  plate 
supported  at  the  corners  for  various  values  of  rw. 

16.  Displacement  of  node  adjacent  to  center  for  the  corner  supported  plated 
and  clamped  plate  for  various  values  of  rw- 

17.  Displacement  w  along  a  radius,  line  B-C  in  Fig.  6,  for  a  thick,  circular 
plate  with  clamped  supports  for  various  values  of  rw. 


18.  Performance  of  the  stabilized  element  for  the  single  element  twist. 

19.  Mesh  for  rhombic  plate  problem;  uniformly  loaded  with  simple  supports 


Zero-Energy  Modes  At  Biquadratic  Plate  Element 
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At  the  present  time,  this  mode  is  valid  for  rectangular  shape  element 
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Boundary  cond. 
u(0,y)  =  0  0<y<l 

u(x,l)  =  0  0<x<0.5 

u(x,0)  =  sin7rx  0<x<0.5 
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Single  Element  Twist 
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